L18: Geometric operations with rasters

Mosaicing, Resampling, Reprojecting

Bogdan G. Popescu

bogdan.popescu@johncabot.edu

John Cabot University

Mosaicing Rasters

In this lecture, we will look at other types of operations that we can perform on rasters:

Mosaicing
resampling
reprojecting

We will continue to use packages like stars and units

Mosaicing: Elevation

We will demonstrate the use the three tools, by applying them to a Digital Elevation Model (DEM).

This is what we try to achieve:

Downloading DEM data

DEM stands for a digital elevation model

There are many data sources

The most helpful website for downloading DEM is: http://www.webgis.com/terr_world.html

Working with Elevation

To save time, you can download the DEM.

You can also download the world borders.

Working with Elevation

Elevation

#Step1: Load the DEM files for every continent
#EUROPE
filepath<-"./data/elevation/w020n90/W020N90.DEM"
x1<-read_stars(filepath)
filepath<-"./data/elevation/w020n40/w020n40.DEM"
x2<-read_stars(filepath)
filepath<-"./data/elevation/e020n40/e020n40.DEM"
x3<-read_stars(filepath)
filepath<-"./data/elevation/e020n90/e020n90.DEM"
x4<-read_stars(filepath)
filepath<-"./data/elevation/w060n40/w060n40.DEM"
x5<-read_stars(filepath)
filepath<-"./data/elevation/w060n90/w060n90.DEM"
x6<-read_stars(filepath)

#ASIA
filepath<-"./data/elevation/e140n40/e140n40.DEM"
x7<-read_stars(filepath)
filepath<-"./data/elevation/e140n90/e140n90.DEM"
x8<-read_stars(filepath)
filepath<-"./data/elevation/e060n40/e060n40.DEM"
x9<-read_stars(filepath)
filepath<-"./data/elevation/e060n90/e060n90.DEM"
x10<-read_stars(filepath)
filepath<-"./data/elevation/e100n40/e100n40.DEM"
x11<-read_stars(filepath)
filepath<-"./data/elevation/e100n90/e100n90.DEM"
x12<-read_stars(filepath)

#AUSTRALIA
filepath<-"./data/elevation/e060s10/e060s10.DEM"
x13<-read_stars(filepath)
filepath<-"./data/elevation/e100s10/e100s10.DEM"
x14<-read_stars(filepath)
filepath<-"./data/elevation/e140s10/e140s10.DEM"
x15<-read_stars(filepath)
filepath<-"./data/elevation/e060s60/e060s60.DEM"
x16<-read_stars(filepath)
filepath<-"./data/elevation/e120s60/e120s60.DEM"
x17<-read_stars(filepath)

#NORTH AMERICA
filepath<-"./data/elevation/w180n90/w180n90.DEM"
x18<-read_stars(filepath)
filepath<-"./data/elevation/w140n90/w140n90.DEM"
x19<-read_stars(filepath)
filepath<-"./data/elevation/w100n90/w100n90.DEM"
x20<-read_stars(filepath)
filepath<-"./data/elevation/w180n40/w180n40.DEM"
x21<-read_stars(filepath)
filepath<-"./data/elevation/w140n40/w140n40.DEM"
x22<-read_stars(filepath)
filepath<-"./data/elevation/w100n40/w100n40.DEM"
x23<-read_stars(filepath)

#SOUTH AMERICA
filepath<-"./data/elevation/w180s10/w180s10.DEM"
x24<-read_stars(filepath)
filepath<-"./data/elevation/w140s10/w140s10.DEM"
x25<-read_stars(filepath)
filepath<-"./data/elevation/w100s10/w100s10.DEM"
x26<-read_stars(filepath)
filepath<-"./data/elevation/w180s60/w180s60.DEM"
x27<-read_stars(filepath)
filepath<-"./data/elevation/w120s60/w120s60.DEM"
x28<-read_stars(filepath)

#SOUTH AFRICA
filepath<-"./data/elevation/w060s10/w060s10.DEM"
x29<-read_stars(filepath)
filepath<-"./data/elevation/w020s10/w020s10.DEM"
x30<-read_stars(filepath)
filepath<-"./data/elevation/e020s10/e020s10.DEM"
x31<-read_stars(filepath)
filepath<-"./data/elevation/w060s60/w060s60.DEM"
x32<-read_stars(filepath)
filepath<-"./data/elevation/w000s60/w000s60.DEM"
x33<-read_stars(filepath)

Elevation

#Step2: Creating a mosaic
star_mosaic_eu <- st_mosaic(x1, x2, x3, x4, x5, x6)
star_mosaic_rus <- st_mosaic(x7, x8, x9, x10, x11, x12)
star_mosaic_aus <- st_mosaic(x13, x14, x15, x16, x17)
star_mosaic_namer <- st_mosaic(x18, x19, x20, x21, x22, x23)
star_mosaic_samer <- st_mosaic(x24, x25, x26, x27, x28)
star_mosaic_safr <- st_mosaic(x29, x30, x31, x32, x33)

#Step3: Creating a second mosaic
all<-st_mosaic(star_mosaic_eu, star_mosaic_rus, star_mosaic_aus, star_mosaic_namer, star_mosaic_samer, star_mosaic_safr)

#Step4: Higher dx/day values means lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4, dy = 0.4))

#Step5: Renaming the raster file
names(elev)<-"elev"

#Step6: Ensuring the the raster has the same crs
borders <- st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet = TRUE)
st_crs(elev)<-st_crs(borders)


fig<-ggplot()+
  geom_stars(data=elev)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM)

If we don’t mosaic the different tiles, we only have one.

This is what the tile for Europe looks like.

Digital Elevation Model (DEM) in Europe

This is how we identify the coordinates of our frame.

#Identifying coordinates of Northern points
#Step1: Selecting the country
norway<-subset(borders, SOVEREIGNT == "Norway")
#Step2: Extracting its coordinates
max_lat_y_eu<-st_bbox(norway)["ymax"]

#Identifying coordinates of Southern points
#Step1: Selecting the country
greece<-subset(borders, SOVEREIGNT == "Greece")
#Step2: Extracting its coordinates
min_lat_y_eu<-st_bbox(greece)["ymin"]

#Identifying coordinates of Eastern points
ukraine<-subset(borders, SOVEREIGNT == "Ukraine")
#Step2: Extracting its coordinates
max_lon_x_eu<-st_bbox(ukraine)["xmax"]

#Identifying coordinates of Western points
#Step1: Selecting the country
portugal<-subset(borders, SOVEREIGNT == "Portugal")
#Step2: Extracting its coordinates
min_lon_x_eu<-st_bbox(portugal)["xmin"]

Digital Elevation Model (DEM) in Europe

#Step1: Higher dx/day values mean lower resolution
all2<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.1, dy = 0.1))
#Step2: Renaming raster values
names(all2)<-"elev"
#Step3: Mapping
fig<-ggplot()+
  geom_stars(data=all2)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.05, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM) in Europe

Digital Elevation Model (DEM) in the US

We can do the same thing for the US

#Step1: Higher dx/day values mean lower resolution
all2<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.1, dy = 0.1))
#Step2: Renaming raster values
names(all2)<-"elev"
#Step3: Identifying the coordinates for the map
min_lon_x_us<-(-130)
max_lon_x_us<-(-57)
min_lat_y_us<-(26)
max_lat_y_us<-(53)
#Step4: Mapping
fig<-ggplot()+
  geom_stars(data=all2)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.05, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_us-3, max_lon_x_us+3), ylim = c(min_lat_y_us-3, max_lat_y_us+3), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM) in the US

We can do the same thing for the US:

Elevation 3D

#Step1: Turning the raster into points
elev2<-st_as_sf(elev, as_points = FALSE, merge = FALSE)
#Step2: Mapping
fig1<-ggplot() +
  geom_sf(data = elev2, aes(fill = elev), color=NA) +
  scale_fill_gradient(low = "chartreuse3", high = "red")+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()
#Step3: Creating the movie
library("rayshader")
library("av")
plot_gg(fig1,multicore=TRUE,width=7,height=3,scale=350)    # Plot_gg de rayshader
render_snapshot(filename = "Elevation")
render_movie("./figures/movie_elev.mp4",frames = 720, fps=30,zoom=0.8,fov = 30)

Elevation 3D

Raster Resampling

Raster resampling is the process of transferring raster values from the original grid to a different grid.

Resampling may be necessary for:

aligning several input rasters that come from different sources to the same grid
reducing the resolution of very detailed rasters, so that they are more convenient to work with

Raster Resampling

In case you have not noticed, we have already resampled the DEM raster.

#Step2: Creating a mosaic
star_mosaic_eu <- st_mosaic(x1, x2, x3, x4, x5, x6)
star_mosaic_rus <- st_mosaic(x7, x8, x9, x10, x11, x12)
star_mosaic_aus <- st_mosaic(x13, x14, x15, x16, x17)
star_mosaic_namer <- st_mosaic(x18, x19, x20, x21, x22, x23)
star_mosaic_samer <- st_mosaic(x24, x25, x26, x27, x28)
star_mosaic_safr <- st_mosaic(x29, x30, x31, x32, x33)

#Step3: Creating a second mosaic
all<-st_mosaic(star_mosaic_eu, star_mosaic_rus, star_mosaic_aus, star_mosaic_namer, star_mosaic_samer, star_mosaic_safr)

#Step4: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4, dy=0.4))

#Step5: Renaming the raster file
names(elev)<-"elev"

#Step6: Ensuring the the raster has the same crs
st_crs(elev)<-st_crs(borders)

fig<-ggplot()+
  geom_stars(data=elev)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Nearest Neighbor Resampling

There are a variety of other re-sampling methods.

In the nearest neighbor sampling, the new raster pixels get the value from the nearest pixel of the original raster.

Some values may be lost.

Nearest Neighbor Resampling

#Step4: Lower dx values. Lower dx values means lower resolution
grid = st_as_stars(st_bbox(all), dx = 0.4, dy = 0.4)
dem1 = st_warp(elev, grid, method = "near")

fig<-ggplot()+
  geom_stars(data=dem1)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig

Bilinear resampling

In bilinear resampling, each new raster cell gets a weighted average of four nearest cells from the input, rather than just one.

Bilinear resampling

Two technical things to note about st_warp when using a method other than method="near":

use_gdal=TRUE needs to be specified, otherwise the method argument is ignored

no_data_value needs to be set to a value outside of the valid range (such as no_data_value=-9999 for a DEM raster)

dem2 = st_warp(elev, grid, method = "bilinear", use_gdal = TRUE, no_data_value = -9999)

Average resampling

Another useful method is the average resampling method, where each new cell gets the weighted average of all overlapping input cells:

Average resampling

dem3 = st_warp(elev, grid, method = "average", use_gdal = TRUE, no_data_value = -9999)

Raster Reprojection

Raster reprojection is more complex than vector layer reprojection.

It requires both transforming pixel coordinates and a resampling step in order to “arrange” the transformed values back into a regular grid.

While other sources recommend st_warp for reprojection, this a lot of pressure on your computer’s memory.

Raster Reprojection

#Step1: Reprojecting the shapefile to EPSG:32636
borders_rep<-st_transform(borders, crs = 32636)
#Step2: Subsetting countries
italy<-subset(borders_rep, SOVEREIGNT=="Italy")
#Step3: Reducing the original resolution: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4, dy=0.4))
#Step4: Turning the raster to a raster object 
dem_rep <- terra::rast(elev)
#Step5: Re-projecting the rasters
dem_rep<-terra::project(dem_rep, terra::crs(borders_rep))
#Step6: Turning the raster back to a stars object
dem_rep<-st_as_stars(dem_rep)

#The following section is to indetify the min max of Italy in ordet to map it.
#Step7: Making a box object around the polygon
bbox <- st_bbox(italy)

#Step8: Extracting min and max latitudes and longitudes
min_lat_y <- bbox["ymin"]
max_lat_y <- bbox["ymax"]
min_lon_x <- bbox["xmin"]
max_lon_x <- bbox["xmax"]

#Step9: Adding an error
error<-10000000
#Step10: Mapping
fig<-ggplot()+
  geom_stars(data=dem_rep)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  #geom_sf(data=italy, linewidth = 0.01, fill = NA, color = "black", alpha=0.5)
  geom_sf(data=italy, linewidth = 0.1, fill = NA)+
  geom_sf(data=borders_rep, linewidth = 0.1, fill = NA)+
  coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Raster Reprojection

fig

Raster Reprojection

Note that the coordinate units in the reprojected raster are no longer degrees, but meters.

We can check that by analyzing the output from st_crs(dem_rep)

st_crs(dem_rep)

Coordinate Reference System:
  User input: WGS 84 / UTM zone 36N 
  wkt:
PROJCRS["WGS 84 / UTM zone 36N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 36N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",33,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 30°E and 36°E, northern hemisphere between equator and 84°N, onshore and offshore. Belarus. Cyprus. Egypt. Ethiopia. Finland. Israel. Jordan. Kenya. Lebanon. Moldova. Norway. Russian Federation. Saudi Arabia. Sudan. Syria. Turkey. Uganda. Ukraine."],
        BBOX[0,30,84,36]],
    ID["EPSG",32636]]

Combining Raster and Vector Layers

Now, we will focus on elements related to:

Cropping and masking a raster according to a vector layer
Switching from vector to raster representation, and vice versa
Calculating a raster of distances to nearest point
Extracting raster values from locations defined by a vector laye

Masking and Cropping Rasters

Masking a raster means turning pixels values outside of an area of interest — defined using a polygonal layer — into NA.

Masking and Cropping Rasters

#Step1: Only choosing Italy
italy<-subset(borders, SOVEREIGNT=="Italy")
#Step2: Reducing the resolution: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.08, dy = 0.08))
#Step3: Cropping
dem_italy<-st_crop(elev, italy)
#Step4: Renaming the raster file
names(dem_italy)<-"elev"

fig<-ggplot()+
  geom_stars(data=dem_italy)+
  scale_fill_gradientn(colours=c("chartreuse3","red"), na.value = NA)+
  geom_sf(data=italy, linewidth = 0.1, fill = NA)+
  theme_bw()
fig

Turning a Raster into a Vector

It is possible to turn rasters into vectors

We have already seen this in fact.

library(sf)
#Step0: Switched off spherical geometry (s2)
sf::sf_use_s2(FALSE)
#Step1: Only choosing Italy
italy<-subset(borders, SOVEREIGNT=="Italy")
#Step2: Ensuring the same CRS for vector and raster
st_crs(italy)<-st_crs(borders)
st_crs(elev)<-st_crs(borders)
#Step3: Cropping raster
dem_italy<-st_crop(elev, italy)
#Step4: Changing the name
names(dem_italy)<-"elev"
#Step5: Creating a grid from the raster
dem_italy2<-st_as_sf(dem_italy, as_points = FALSE)
#Step6: Changing the name of the raster
names(dem_italy2)[1]<-"elev"

Turning a Raster into a Vector

This fundamentally creates a grid

#Step7: Making a box object around the polygon
bbox <- st_bbox(italy)

#Step8: Extracting min and max latitudes and longitudes
min_lat_y <- bbox["ymin"]
max_lat_y <- bbox["ymax"]
min_lon_x <- bbox["xmin"]
max_lon_x <- bbox["xmax"]

error<-1
fig1<-ggplot()+
  geom_sf(data=dem_italy2, aes(fill=elev), linewidth=0)+
  geom_sf(data=italy, fill=NA)+
  ggtitle("Shapefile")+
  coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

fig2<-ggplot()+
  geom_stars(data=dem_italy)+
  geom_sf(data=italy, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Raster")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Turning a Raster into a Vector

This fundamentally creates a grid

library(ggpubr)
ggarrange(fig1, fig2)

Turning a Raster into a Vector

Extracting Raster Values

It may be necessary to “extract” raster values according to:

Extracting Raster Values

It may be necessary to “extract” raster values according to:

overlapping points

Extracting Raster Values

It may be necessary to “extract” raster values according to:

overlapping line

Extracting Raster Values

It may be necessary to “extract” raster values according to:

overlapping polygon

Extracting Raster Values

We can aggregate and take the average of the raster pixels at the polygon level

You can download the shapefiles for Italy.

library(sf)
library(dplyr)
#Step1: Reducing the resolution: Higher dx/day values mean lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.2, dy = 0.2))
#Step2: Renaming the raster
names(elev)<-"elev"
#Step3: Reading the shapefile
ita2 <- st_read(dsn="./data/gadm41_ITA_shp/gadm41_ITA_2.shp", quiet = TRUE)
#Step4: Keeping only the relevant variable
ita2<-subset(ita2, select = c(NAME_2))
#Step5: Simplify polygons
ita2<-st_simplify(ita2,  dTolerance = 0.005)
#Step6: Making sure that the raster and shapefile have the same crs
st_crs(ita2)<-st_crs(elev)
#Step7: Cropping the raster to the shapfile
dem_italy<-st_crop(elev, ita2)
#Step8: Calculate the mean by polygon
elev_by_county <- raster::aggregate(dem_italy, ita2, mean, na.rm=TRUE)
#Step9: Adding the raster by polygon
poly_elev <- mutate(ita2, avg_elev = elev_by_county$elev)

Extracting Raster Values

We can for example aggregate and take the average of all rasters at the polygon level

Loading counties in Italy

fig1<-ggplot()+
  geom_sf(data=ita2, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Polygons")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

fig2<-ggplot()+
  geom_stars(data=dem_italy)+
  geom_sf(data=ita2, fill=NA)+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Raster")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

fig3<-ggplot()+
  geom_sf(data = poly_elev, mapping = aes(fill = avg_elev))+
  scale_fill_continuous(na.value=NA)+
  ggtitle("Zonal Stats")+
    coord_sf(xlim = c(min_lon_x - error, max_lon_x + error), 
           ylim = c(min_lat_y - error, max_lat_y + error))

Extracting Raster Values

Rasterizing Polygons

We might be interested in turning the shapefiles that we just create into a raster.

We can achieve that by using st_rasterize

poly_elev_raster<-st_rasterize(poly_elev)

ggplot()+
  geom_stars(data=poly_elev_raster)+
  scale_fill_continuous(na.value=NA)

Conclusion

We have learned a variety of additional raster tools:

working with DEM (Digital Elevation Models)
mosaicing, resampling, reprojecting rasters
masking and cropping rasters
turning a raster into a vector and a vector into a raster